Simulation of the Sedimentation of a Falling Oblate 
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We present a numerical investigation of the dynamics of one falling oblate ellipsoid particle in a 
viscous fluid, in three dimensions, using a constrained- force technique [19], [20] and [17]. We study 
the dynamical behavior of the oblate for a typical downward motion and obtain the trajectory, 
velocity, and orientation of the particle. We analyze the dynamics of the oblate generated when the 
height of the container, the aspect-ratio, and the dynamical viscosity are changed. Three types of 
falling motions are established: steady-falling, periodic oscillations and chaotic oscillations. In the 
periodic regime we find a behavior similar to the case of falling flat strips reported in ref. [13]. In the 
chaotic regime the trajectory of the oblate is characterized by a high sensitivity to tiny variations 
in the initial orientation. The Lyapunov exponent is A = 0.052 ± 0.005. A phase space comparing 
to the results of ref [12], is shown. 



I. INTRODUCTION 

The way in which objects fall to the ground has been 
studied since antiquity. Objects were thought to return 
to "their natural" places by the ancient Greeks. Newton 
showed that the bodies fall on earth driven by a constant 
acceleration. But despite gravity's undeniable attraction, 
not all falling objects travel downwards on straight tra- 
jectories. The tree leaves flutter to the ground in the 
autumn, exhibiting a complex motion and refusing to 
follow the shortest path. 

A deep understanding of the motion of falling objects 
in a fluid is of great technical importance, and has been 
investigated in a variety of contexts, including meteorol- 
ogy [1] , aircraft stability [2] , power generation [3] , chem- 
ical engineering [4], etc. Also Newton observed the com- 
plex motion of objects falling in both air and water [6]. 
This phenomenon was also studied by Maxwell, who dis- 
cussed the motion of a falling paper strip [7] . 

In the nineties, Aref and Jones [10] found through nu- 
merical solutions of Kirchhoff's equations that the tra- 
jectory for an object moving through an incompress- 
ible inviscid and irrotational fluid, is chaotic. Tanabc 
and Kaneko [9], using a phenomenological model for the 
falling of a one-dimensional (ID) piece of paper, including 
lift and dynamical viscosity, but neglecting the inertia of 
the fluid, describe five falling regimes. Two of them are 
chaotic. In 1997 Stuart B. Field et al. [12], investigated 
experimentally the behavior of falling disks in a fluid, 
and identified different dynamical regimes as function of 
the moment of inertia and the Reynold's number. They 
obtained experimental evidence for chaotic intermittency 
[11]. In 1998 Andrew Belmonte, et al. [13], in an exper- 
iment with thin flat strips falling through a fluid, ob- 
served only two motions: side to side oscillation (flutter) 
and end-over-end rotation (tumbling). They proposed a 
phenomenological model including inertial drag and lift 
which reproduces this motion, and yields the Froude sim- 
ilarity, which describes the transition from flutter to tum- 



ble regime. Mahadevan et al [23], in 1999 made an ex- 
periment of dropping horizontal cards of thickness d and 
width w, showing that the tumbling frecuency f2 scales 
as 17 rs c? 1 / 2 w~ 1 , consistent with a dimensional argument 
that balances the drag against gravity. 

Given the difficulties to study this problem theoreti- 
cally and experimentally, we took a computational ap- 
proach simulating the falling of one oblate ellipsoid in a 
viscous fluid in a three dimensional container. An oblate 
is an ellipsoid for which the two largest principal radia 
are equal. We organize the paper in the following man- 
ner. In section 2 we give an review over the model that 
we use. In Sec. 3A we describe the main features of the 
falling oblate. In Sec. 3B-D we present the results for the 
change in the initial height, the oblate's aspect-ratio and 
the dynamical viscosity, for the steady-falling regime. In 
sec. 3E we show the periodic behavior and compare to 
the results of reference [13]. In Sec. 3F the sensitivity to 
tiny variations in the oblate's initial orientation is pre- 
sented in the chaotic regime. In sec. 3G the parameter 
phase space is sketched and compared to the results of 
reference [12]. Section 4 summarizes present results and 
discusses possible further applications. 



II. MODEL 



The general idea, proposed by Fogelson and Peskin 
[21], is to work with a simple grid for the resolution of 
the fluid motion at all times and represent the particles 
not as boundary conditions to the fluid, but by a volume 
force term or Lagrange multipliers in the Navier Stokes 
equations. 

This technique was developed in the work of Kuusela, 
et al [17], Wachmann, et al [18] and Hoefler et al. [19], 
[20], and employs a numerical solver for the dynamical 
simulation of three-dimensional rigid particles in a New- 
tonian fluid, bounded by a rectangular container. 
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The motion of the fluid is described by the dimension- 
less Navier-Stokes equations: 

^ + (v-V)v = -Vp+±-V 2 v + f (1) 

Here p and v are the pressure and the velocity of the 
fluid, respectively, and / is an external force. 
The Reynolds number Re is defined as 



Re := 



UDp f 



(2) 



where U is the mean vertical oblate velocity, D a char- 
acteristic length that in our case is the largest oblate's 
diameter, pf the density and v the dynamical viscosity 
of the fluid. 

For an incompressible fluid, the continuity equation: 



dpj_ 
dt 



reduces to 



+ V • (p f v) = 







(3) 



(4) 



Equation (1) is discretizcd on a regular, marker-and- 
cell mesh to second order precision in space. For the 
time stepping, we employ an operator-splitting-techniquc 
which is explicit and accurate to first order. The suspend- 
ing fluid is subjected to no-slip boundary conditions at 
the surface of the suspended particles. More details of 
the solution procedure are presented in [17], [19], [20]. 

The geometry of the oblate ellipsoid is characterized 
by Ar its aspect-ratio defined as the ratio of the smallest 
radius over the largest one. 

An oblate is represented by a rigid template connected 
to fluid tracer particles, which are moving on the tra- 
jectories of the adjacent fluid. The connection is made 
by using the body force term, in the Navier-Stokes equa- 
tions, as constraints on the fluid such to describe the 
oblate. 

The force density / c , is chosen elastic with a spring 
constant that, and guarantees that the elongation re- 
mains small against the grid spacing at all times [18], 
and it is zero in the exterior of the region outside the 
oblate. We can define this force density f c as: 



(5) 



where the displacement field of the separation be- 

tween the markers i and their corresponding reference 
point j. The stiffness constant k, must be chosen large 
enough so that | e(x~ij) |<C h, (h size grid), holds for all 
iterations. 

In general the displacement field e{xlj) is defined as: 



£j {xij ) 



(6) 



The vector x7j™ is the position of a fluid tracer, whose 
motion is determined by the fluid local velocity, i.e., 



= u{x? j ) 



(7) 



The xlf are the reference points associated to a tem- 



plate having the shape of the physical particle: 



(8) 



Here X{ is the center of mass of the template, Oi(t) is 
the rotation matrix that describes the present orienta- 
tion of the oblate and r?j denote the initial position of 
the reference points with respect to the center of mass. 
For the quaternion formulation of the rotation, we use 
the techniques described in ref. [15]. 

A velocity- Verlet integrator [16] serves to integrate 
the equations of motion for the translation and a Gear- 
predictor integrator [5] for the rotation on the template: 



F = -Mgj+ Pf Vgj+J2ff 



(9) 



where j is the unit vector along the vertical. 



T — ^ ^ ( x i x cm) x fl 



(10) 



with respect to the template's center of mass f cm . 
The equations of motion of the particle template are: 



(11) 



and 



n 



T 

7 



(12) 



where M is the mass of the template particle; U and f2 
are the linear and angular velocities of the template par- 
ticle, respectively; I is the moment of inertia; and T is 
the torque, [17], [19]. 

The boundary conditions at the container wall are zero 
for the normal velocity component of the fluid and no- 
slip condition for the tangential component, because the 
walls are assumed impenetrable, [19], [18]. The interac- 
tion between the oblate and the walls is defined through 
a contact force, [17], where the walls are treated as a 
particle with infinite mass and infinite radius. 
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III. RESULTS 



Change of Vertical Velocity 




FIG. 1. (a) Steady- falling, (b) Periodic-oscillation, (c) 
Chaotic motion. 



FIG. 2. (top) Decreasing amplitude of the vertical velocity, 
(bottom) Comparison between the vertical velocity and the 
spatial trajectory at the same height. Initial conditions of the 
system. 8 a = 26.6°, Ar = 0.25, ^ = 0.033. Falling initial 
height h a = 76cm for the case of steady-falling. 



For all our simulations, and in order to reduce the pa- 
rameter space we use p f i uid = 1^ and p ob iate = 3.5^. 

We found three basically different motions in our sim- 
ulations: steady-falling, side-to-side periodic-oscillation 
known as 'flutter' ref. [13], and a chaotic motion as shown 
in fig. 2. The above phenomenology can be compared 
to the work of ref. [12], for the case of dropping disks. 
In general, the trajectories depend strongly on the ini- 
tial conditions and the properties of the system (oblate's 
orientation o , dynamical viscosity /x and the oblate's 
aspect-ratio Ar, etc). We don't find the tumbling mo- 
tion described in the above references. 



A. Phenomenology of the Steady-Falling Oblate 



For the velocity of the center of mass, the vertical com- 
ponent decreases when the oblate approaches the con- 
tainer bottom and shows a damped wavering, with an 
amplitude that decreases with time (fig. 3 top). 

The vertical trajectory is composed of succesive turn- 
ing points, that correspond to the points where the tra- 
jectory changes the sign of the rate of change of the verti- 
cal velocity component marked in fig. 3 (bottom) by the 
horizontal lines L2, ... L5. We also see, that from L\ 
to Lg the amplitude for the trajectory and the velocity 
decrease. 
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Horizontal Velocity & T 
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FIG. 3. Horizontal velocity. Initial conditions of the sys- 
tem. 6 = 26.6°, Ar = 0.25, fi = 0.033. Falling initial height 
(b)h — 76cm. 

The horizontal components of the velocity obey a dif- 
ferent behavior as the vertical ones, and in general they 
have non-regular oscillations as seen in fig. 4. 



The oblate orientation is described through the three 
rotational degrees of freedom, around the center of mass. 
We present the time evolution of the angle between the 
oblate's normal and the vertical direction that we will 
call vertical orientation fig. 5 (top). 6 = implies that 
the oblate's principal axis will be parallel to the con- 
tainer's horizontal fig. 5 (bottom). At the beginning of 
the movement, there is a larger angular change of the 
oblate's normal AO fig. 5(bottom), which is character- 
ized by a larger peak-to-peak amplitude G p _ p . In fig. 5 
the definitions we illustrate the definition of the peak-to- 
peak amplitude as the distance between succesive turn- 
ing points, which decreases while the oblate sinks. In the 
steady-falling regime, both the peak-to-peak amplitude 
Op_ p and the change A9 = 0/ — ©j , become smaller 
from the top towards the container's bottom. The oblate 
tends to align its major axis along the vertical [14], pre- 
senting the lowest resistance to its descent in the fluid, 
and acquiring a limit vertical velocity fig. 3 (top). 



Change of Vertical Orientation 




15 Decreasing peak-to-peak amplitude 

°0 5 10 15 20 25 30 35 40 
Time t 

Angular Variation 

85 1 ■ ■ ■ ■ ■ ■ 1 




Position in X 



FIG. 4. (Top) Decreasing peak-to-peak amplitude in the 
angular oscillation. (Bottom) Angular change AO along 
the vertical trajectory. Initial conditions of the system. 
9 = 26.6°, Ar = 0.25, fi = 0.033. Falling initial height 
(b)h — 76cm. 




FIG. 5. The right picture shows the vortex structure of the 
vertical and horizontal components of the fluid velocity field 
(it, v) generated by the falling oblate, with diameter of 3.2 
cm in a container of 10 x 30 x 30 cm and Reynolds number 
of Re = 128, aspect-ratio Ar = 0.5. The left picture shows 
shedding vortices reported by Belmonte et al, ref.[12], for a 
falling strip. 



4 



The oblate generates shedding vortices in the fluid 
along its falling trajectory, an example is shown in figure 
6 (right), that shows the velocity fluid field around the 
oblate, and the vortex is localized just in the top region 
above the oblate, and where the oblate has experienced 
the larger angular change AO as shown in fig. 5 (top). 
The Reynolds number calculated from the oblate's di- 
ameter and terminal velocity is Re = 128. We point out 
that the vortex structure is obtained also in the work of 
Belmonte et al, [13] where a shedding vortex created by 
the zigzag motion of the falling strip is seen in the left 
figure 6. 



Steady-Falling Oblate: Change in the Initial 
Height. 



Vertical Orientation vs. Time 




15 20 25 
Time T 



Oblate Vertical Trajectory 
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Horizontal Position 



Vertical Velocity vs. Time 
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FIG. 6. Initial conditions of the system. O — 26.6°, 
Ar = 0.25, /i = 0.033. Each Trajectory has different initial 
height (a)h = 96cm, (b)h = 76cm, (c)h = 56cm. I) The 
spatial trajectory in the vertical plane. II) Vertical velocity 
vs. time. Ill) Vertical orientation vs. time. 



We choose three different initial heights fig. 71, (fixing 
the rest of parameters). If we superpose the trajectories 
in fig. 71, it can be shown that there is no variation in the 
wavelength or in the peak-to-peak amplitude. For all the 
heights, the trajectories generated are in good agreement 
with a damped harmonic oscillation. 



For the range of falling heights used, the oblate at- 
tains the same terminal vertical velocity fig. 711, (3^), 
its magnitude is independent on the falling height. The 
vertical velocity finally, reaches a stable state (uniform 
and linear motion), after the same time (~ 20s fig 711). 
We also see that the vertical velocity suddenly becomes 
zero when the oblate touches the bottom. 



We see in fig. 7 Ill(c-a), an increase in the final angle 
with respect to the initial height. For h Q — 56 the smaller 
height, we obtain ~ 85° and for larger one h = 96, we 
have ~ 95°. For the larger height the oblate is still in 
a transitory state before it arrives to its final angle fig. 
7III. 
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C. Steady-falling oblate: Change in the dynamical 
viscosity. 



Oblate Vertical Trajectory 

90, ■ ■ ■ ■ 




TimeT 



FIG. 7. Initial conditions in the system. o = 26.6*, 
h — 80, Ar = 0.25. Each trajectory has a different dy- 
namical viscosity (a)fi = 0.025, (b)fi = 0.033, (c)/tt = 0.100. 
I) Proyection of the trajectory in the vertical plane. II) Ver- 
tical (top) and horizontal (bottom) velocities vs. time. Ill) 
Vertical orientation vs. time . 



The oblate starts its swinging motion with a given ini- 
tial orientation 8 = 26.6°. It glides downwards and to 
the side acquiring some amplitude, while the dynamical 
viscosity n acts reducing the subsequent amplitudes of 
oscillation (fig. 81) . As we increase the dynamical viscos- 
ity from (o)/i = 0.025 to (c)/i = 0.100 the attenuation 
in the oscillatory trajectory becomes stronger. In fig. 81 
a, there is a weak damping producing a long oscillatory 
behavior and in fig. 81 c, the trajectory is quickly atten- 
uated in the first half of the falling height, and hereafter 
it follows a vertical trajectory Similar decreasing oscil- 
lations, are observed when one small sphere suspended 
from a fixed point by a string or a rod, oscillates against 
the air and generates a damped harmonic oscillator. 



The behavior of the vertical velocity is shown in fig. 
811. There is a clear damping in the peak-to-peak velocity 
amplitude, in all the cases. The attenuation in the ver- 
tical velocity amplitude and the time between two con- 
secutive turning points are not very different from each 
other, as the fluid dynamical viscosity is changed. For 
the three viscosities the oblate adquires approximately 
the same final mean velocity v y — 3—. 



A larger variation in the velocity is observed in its hor- 
izontal component fig. 811. There is a strong attenuation 
in the curve for larger viscosities fig. 8IIc. One explana- 
tion is that the interaction between the walls and the 
oblate is smaller when the viscosity is decreased [22]. 
This is also supported by the larger angular variation 
A 9 in fig. 8IIc, for smaller viscosities. 



In fig. 8 III we see that the first peak-to-peak ampli- 
tude are nearly the same for all considered viscosities. 
For fj, — 0.025 the subsequent peak-to-peak amplitudes 
attain a constant value of 6 p _ p ^15° while for higher 
viscosities periodic oscillations arc replaced by erratic 
motion of low amplitude. 
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D. Steady-Falling Oblate: Change in the 
Aspect-Ratio of the Oblate. 



Increment in the Aspect-Ratio 




23456789 10 
Horizontal Position 

Vertical Velocity vs. Time 




FIG. 8. Initial conditions of the system. O = 26.6°, 
h = 96, L x x L y x L z = 10 x 100 x 10, n = 0.033. Each trajec- 
tory has different aspect-ratio (o)Ar = 0,29, (fe)Ar = 0,22, 
(c)Ar = 0, 18. I) Trajectory for the vertical plane. II) Vertical 
velocity vs. time. Ill) Vertical orientation vs. time. 



When the oblate's falling motion begins, the oblate 
gains a larger oscillation amplitude in its oscillatory 
falling trajectory. And for all the three trajectories pre- 
sented in fig. 9, this first amplitude decreases A Q = 
3.0,2.8,1.8cm as the oblate's aspect-ratio is increased 
Ar = 0.18,0.22,0.29 respectively. This first large oscil- 
lation is a common characteristic for these trajectories. 
As the aspect-ratio is incremented, the number of cy- 
cles and the amplitude of the trayectories change. For 
Ar = 0.22 the trajectory presents a well defined steady- 
falling behavior, fig. 91(b), but, if the aspect-ratio in- 
creases, (Ar = 0,29), the trajectory varies, the peak-to- 
peak amplitude is quickly damped in the first half of the 
vertical trajectory, with the interesting observation that 
in the second half the trajectory doesn't have a steady- 
falling behavior fig. 9 la. When Ar = 0.18, the smaller 
aspect-ratio, the trajectory has a oscillatory behavior fig. 
91(c), with a constant peak-to-peak amplitude of 3cm. 



Since the minor oblate radius is fixed in our sim- 
ulations, when the aspect-ratio is increased Ar = 
0.18,0.22,0.29, the oblate's area gets smaller, and the 
final vertical velocity increases Vy — 3.5,3.7,3.9^, re- 
spectively, fig. 911. The final vertical velocity decreases 
with the decrement in the aspect-ratio, since the smaller 
aspect-ratio presents more area against the fluid. As the 
aspect-ratio is increased the peak-to-peak amplitude in 
the vertical velocity diminishes. For the smaller aspect- 
ratio Ar = 0.18, we have the larger amplitude l^p, and 
as the aspect-ratio is increased Ar = 0.22 and Ar = 0.29 
the peak-to-peak amplitude tends to be much smaller. 



15 20 
Time T 



The peak-to-peak amplitude for the vertical orienta- 
tion Qp-p increases when the aspect-ratio decreases or 
the oblate's area increases. For Ar = 0.29 the peak- 
to-peak amplitude is ® p - p = 15°, (fig. 9a), and much 
smaller compared to Q p - P = 70°, (fig. 8c), for Ar = 0.18. 
In all cases the oblate at the end orients vertically (fig. 9 
III). 
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E. Periodic Behavior of a Falling Oblate. 



Oblate Trajectory 
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Vertical Orientation & Time 
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H=0.025 



H=0.100 




15 20 25 
Time 

FIG. 9. Trajectories generated for /Hi = 0.100, /t*2 = 0.025. 
(I) Trajectory in vertical and horizontal position. (II) Verti- 
cal orientation O vs time. (Ill) Vertical velocity v y vs time. 
The initial conditions are h a = 96, Ar = 0.133, 6 = 63.4°. 

We have found periodic behavior for smaller dynami- 
cal viscosity (Re = 480) and smaller aspect-ratio (Ar = 
0.133). The dynamics of the falling oblate is governed 
by incrtial effects. In figure 10 (I), wc show the tran- 



sition from a quasi-periodic, or a long steady-falling 
trajectory (/i! = 0.100), to a periodic behavior fig. 10 
(I,/i2 = 0.025), when the dynamical viscosity is varied 
from hi = 0.1 to [ii = 0.025. The trajectory presented 
in fig. 101, with dynamical viscosity \±i — 0.025 has a 
wave length of 20 cm. 

The vertical velocity shown in fig. 10 (II), has the same 
transition from a long steady-falling regime with a final 
average velocity of 3.0^ to the periodic regime where 
the velocity has a oscillation period of 3.3 s. 

The vertical orientation presented in fig. 10 (III) has 
also the same transition from a long steady-falling regime 
to periodic behavior with a period of 6.6 s, and the angu- 
lar values oscillate around 9 = 90° with angular peak- 



to-peak amplitude 6 



p-p 



60°. 
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Vertical Orientation & Time 




5 10 15 20 25 30 35 40 
TimeT 



FIG. 10. I) Trajectories for three initial orientations, (a) 
Go = 26° (b) Go = 63° (c) 9 = 90° . II) The corresponding 
vertical velocities. Ill) The vertical orientations. Initial con- 
ditions h a = 100, L x x L z = 10 x 10, Ar = 0.133, Re = 435. 
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a) Simulation 
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We perform three simulations in the periodic regime 
with very different initial orientation angles and the cor- 
responding trajectories are shown in fig. 111. In the case 
of 0o = 26° the peak-to-peak amplitude is 2.3 cm and 
for (Oo = 90°) it is 0.4cm, and the oscillatory behavior 
is observed for the three cases. The peak-to-peak ampli- 
tude of the oscillation in the trajectory fig. 101, decreases 
as O is increased. In the case of the vertical velocity 
and orientation fig. 11II-III, the initial orientation angle 
also plays the same role, reducing the peak-to-peak am- 
plitude of the curves, and for 6 = 90° the amplitude of 
oscillation is the smallest of the three. 

The final vertical velocity and orientation for the three 
trajectories are 3.0^ and 85° respectively. We can say 
that the average final values for the vertical velocity and 
orientation are not modified by the variation of the initial 
orientation 6 - 

We have found the largest peak-to-peak amplitude 
l.O^p for the vertical velocity oscillation for an initial 
orientation of 6 = 26° and it becomes smaller as the 
oblate's initial orientation tends to 9o = 90°. For the 
peak-to-peak amplitude of the vertical orientation the 
oblate shows a similar behavior: we have the larger value 
for the amplitude 6 



p-p 



70°, and it is obtained for an 
initial orientation of 6o = 26°. The smallest peak-to- 
peak amplitude 6 P _ P ~ 4° is obtained for 6q = 90°. 
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F. Sensitivity to a Change in the Initial Orientation 
In the chaotic Regime. 
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FIG. 12. Initial conditions h„ = 96, Ar = 0.25, n = 0.033 
and tiny variations of the initial orientation (a)0 — 26.6°, 
(b)0 o = 26.6001°, (c)0„ = 26.6000001°. 



FIG. 11. Comparison with the results of Belmonte et al. 
ref[12] for a falling strip in the periodic regime, for the (a) 
Vertical orientation Q vs time, (b) Vertical velocity V y vs O. 
(c) Horizontal velocity V x vs O. The initial conditions are 
h = 96, Ar = 0.133, m = 0.025, = 63.4°. 

In figure 12a (simulation) we show the time dependence 
of the vertical orientation with ^ 2 = 0.025. The value of 
the angular peak-to-peak amplitude is 9 = 60°. The 
vertical velocity fig. 12b (simulation) reaches its max- 
imal value 3.6^ as approaches Q ma x, presenting a 
minimal drag. 

The smaller vertical velocity V y = 2.5^p in the oscil- 
lation is reached at <d m in ~ 105°. The butterfly shape 
of fig. 12b was also measured in the experimental work 
of Belmonte et al. [13] (fig. 12b), exhibiting a vertical 
orientation 6 that oscillates at double the period of V y . 

The horizontal velocity oscillates around zero with the 
same period as 0, presenting its maximum value at 
Vx max = 1.5^p and the minimum at Vx m i n = —1.5^ 
at 6 ~ 90° as seen in fig. 12c. When the horizontal 
velocity is zero the oblate takes its maximum (120°) and 
minimum (60°) values in corresponding to a non-zero 
value of the vertical velocity V y = 3.2^p. 



We can ask what is the sensitivity of the oblate to 
tiny changes in the initial orientation. Therefore we have 
simulated three trajectories shown in fig. 13, which have 
slightly different initial orientation. A tiny variation in 
the relative orientation (A0 O = 10~ 3 ) produces, how- 
ever, a significant variation in the shape of all the curves. 
These can be appreciated in the lower part of the trajec- 
tories. 



Trajectory Sensitivity to the Initial Orientation 
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FIG. 13. Initial conditions h a = 166, Ar = 0.25, fj, = 0.033 
and tiny variations of the initial orientation (a)8o = 45.384°, 
{b)0„ = 45.033°, (c)0„ = 44.981°, (d)0 o = 44.976°. 
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In order to get better sensitivity, we have incremented 
the falling height to h Q — 166cm. The resulting trajec- 
tories for four slightly different initial orientations in the 
vertical plane are presented in fig. 14. In this regime 
the system presents a high sensitivity to the initial ori- 
entation condition. For the four trajectories the relative 
angular variation is A0 o = 1CT 3 . 

Autocorrelation Function 
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FIG. 14. Detection of chaos. I) Autocorrelation function 
for the time series of x(t) for the trajectory of fig. 14a. II) 
Power Spectra of fig. 14a III) Poincare section (p x ,x) for the 
trajectory of fig. 14a. 

Due to this dependence on small changes in the initial 
orientation, we proceed to use as a tool of diagnosis, the 
Fourier power spectrum of the time series of the horizon- 
tal coordinate x(t), x(t+St),x(t+2*St)..., and in our case 



5t = 0.053566. A broad spectrum of frecuencies appears, 
as shown in fig. 1511, indicating chaotic motion. 

The autocorrelation function, for the same time series 
(see fig. 151), does not fall quickly to zero, it decreases 
linearly with time. The points are not independent of 
each other and a self similarity is present in the data. 

In the figures 14III, we present a slices or Poincare sec- 
tions (p x , x), corresponding to the trajectories in fig. 14a, 
and which are quite irregular. 

The orbits are quasi-periodic in the sense that they 
pass repeatedly and irregularly through the whole do- 
main without ever closing on themselves, and without 
any particular time period associated with succesive pas- 
sages. 

The sensitivity to initial conditions is clear in these 
four figures. A small change in the initial orientation 
results in large changes in position and velocity. 

Increment of the Distance in Time 
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FIG. 15. Increasing logarithmic behavior for the separation 
distance between the trajectories (a) and (c) in fig. 14, that 
slightly differ in the initial orientation angle by AO = 0.403°. 



We can now investigate quantitatively this sensitiv- 
ity by studying the increment in the Euclidean distance 
^piP2 — V ( x i ~ x 2) 2 + (yi — J/2) 2 , between the curves 
presented in fig. 14 (a) and (c). Fig. 16, shows that 
the distance between nearby points has an overall expo- 
nential time dependence d(t) ~ exp(Ai) and the fit gives 
an estimate for the Lyapunov exponent A = 0.052±0.005. 
The positivity of the Lyapunov exponent is a clear indi- 
cation for chaos. 



G. Parameter Phase Diagram. 

We explore the phase space in the dimensionless mo- 
ment of inertia /*, which is the ratio of the moment of 
inertia of an oblate around its principal axis to the mo- 
ment of inertia of a sphere of liquid with the same di- 
ameter and the Reynolds number Re. We do a similar 
analysis for our results as in the work of Field, et al. [12]. 
It is important to remark that the mentioned experiment 
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was for a falling disk, with small aspect-ratio, and we ex- 
pect that the dynamics of the system will be close to that 
of an oblate ellipsoid. 

The definitions of the dimensionless variables for our 
system are: 



I* = 



'-oblate 



5 r m Poblate 5 poblate 



-'-sphere 4 r M P fluid 4 p fluid 



Ar 



Re = 



U(2r M )Pfiuido 
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Phase Diagram for the Falling Oblate Ellipsoid 
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FIG. 16. The top picture shows the phase diagram of falling 
disks reported in ref. [11]. In the bottom plot we present the 
regimes of the phase space for the falling oblate obtained in 
our simulations. 



Figure 17 bottom, shows our results. At low values 
of I* and small Reynolds number (high dynamical vis- 
cosity), the left-down corner of the diagram, the motion 
is overdamped and the oblate drops to the bottom con- 
tainer without any oscillation. If the Reynolds number 
increases (Re > 100), fixing the moment of inertia, the 
trajectory is composed of successive oscillations that will 
decrease in amplitude until the oblate finally comes to 
stop at the lowest part of the container part. For all 
these cases this type of motion was called in our results 
a steady-falling, and was studied in section IIIB-B,C,D. 

For small values of I* ~ Ar <C 1, we have a flattened 
ellipsoid, and Reynolds number (Re > 400), the trajec- 
tory, velocity and orientation are characterized by oscil- 
lations that repeat at equal intervals of time and space. 
This case was called in our simulations the oscillatory 
regime, and was studied in section IIIB-E. 

As we increase /*, the object will become a sphere 
slightly flattened in the poles, and its dynamic becomes 
sensitive to smaller variations in the initial orientation, 
exihibiting a chaotic trajectory, which is explained in sec- 
tion IIIF. 

If we compare our diagram with the experimental re- 
sults obtained by S. Field et al, ref [12] (fig. 17 top), we 
can see that in both pictures, the steady-falling region 
is located in the lower- left part, the oscillatory region 
is in the lower- right, and the chaotic regime is located 
along the top of the picture. We can say that the two 
diagrams are similar, but, with the difference, that the 
tumbling regime in the Field's diagram is not present in 
our results because by the comparison of the two dia- 
grams in fig. 17. The tumbling regime would be found 
for larger dimensionsless moment of inertia, in our case 
smaller aspect-ratio, demanding a larger three dimen- 
sional container and impliying a very expensive compu- 
tational study, than in the other regimes. 

If we take a different initial orientation angle, the new 
phase diagram exhibits the same dynamical behavior of 
figure 17 bottom. The coexistence of the dynamical 
phases, explained above, is independent of the initial ori- 
entation of the oblate. 



IV. CONCLUSIONS AND OUTLOOK 



The motion of a single oblate settling in a fluid in 
a three dimensional container has been studied. We 
found three basic regimes for the dynamics of the system 
(steady- falling, oscillatory, and chaotic). The steady- 
falling and the periodic motion exhibit a similar physical 
behavior as observed for flattened bodies [12], [13]. With 



12 



the exception that the tumbling motion is missing in our 
simulations. 

We have characterized the dynamics of the steady- 
falling regime when the dynamical viscosity dropping 
height, and oblate's aspect-ratio are changed. Several 
conclusions can be drawn from this part of the work. 

(a) The spatial trajectories (x, y) are composed of os- 
cillations that correspond to a damped harmonic motion. 
This regime is present for small values of I* « 0.5 — 1, 
Re «s 100 and is shown in fig. 6-7. There is no variation 
in the trajectories when we increase the initial height. 
When the aspect-ratio is varied, the trajectories change 
very much, fig. 8. The aspect-ratio dominates strongly 
the type of trajectory that is present in the system. 

(b) The final vertical velocity V y does not depend on 
the initial falling height and the dynamical viscosity. 

(c) The vertical orientation 9 of the oblate, under- 
goes a rotational motion till its mayor axis is aligned 
with the direction of gravity. This tendency of finding 
the minimal resistance against the fluid, is present for all 
Reynolds numbers in the range Re w 30 — 100 used in 
our simulations. 

The periodic behavior in our simulations is found for 
(Re ~ 500), and small (I* < 0.5). The most important 
characteristic on this regime is that the vertical orien- 
tation O, oscillates at the double of the period of the 
vertical velocity V y , and at the same period of the hori- 
zontal velocity V x . This periodic motion is also present 
in the work of Bclmonte et al. [13]. In this regime the 
initial orientation determines the value of the oscillation 
amplitude in the spatial trajectory (x, y), velocity V y and 
orientation 9. For 9 Q = 90° the amplitudes of the above 
quantities approach a smaller value. 

The chaotic behavior is present for I* > 1 and in the 
entire range of Reynolds numbers used in the simula- 
tion. The separation between the spatial trajectories of 
the falling oblate will diverge for small variations in the 
initial orientation 9 Q , and grows exponentially in time. 
The value found for the Lyapunov exponent is around 
A = 0.052 ± 0.005. 

The construction of the phase diagram shows three 
well-defined dynamical regions as in the case of ref. [12]. 
But with the difference that the chaotic behavior in the 
above reference is associated with a transition to chaos 
through intermittency for which we have no indication in 
our simulations. The phase diagram is independent on 
the initial orientation. 

More work to better understand the role of the fluid 
pressure and velocity fields as well as a more systematic 
study of the phase transition in the phase diagram seem 
necessary 
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